
rm(list=ls())

## requires JAGS version 2.1.0, rjags version 2.1.0-2, and R2jags version 0.02-07
## it is incompatible with newer versions of JAGS e.g. 2.2 and 3.x

library(R2jags)
library(foreign)

## dummy variables for missing data
dummies <- TRUE

## choice of violence experience variable
violence.exp <- "area"

load("endorseSetup.RData")

## I'm going to run this three times for APSR
i <- 1

##########################
##  The Village Model  ##
##########################
data <- list ("N", "J", "K", "Y", "endorser", "village", "district",
              "province", "n.village", "n.district", "n.province", "indivcov",
              "villcov", "distcov")
inits <- function() {list(
                          alpha0 = matrix(seq(-2,2,length.out=(4*J)),J,4), 
                          beta0 = runif(J),
                          x = rnorm(N),
                          s = endorser - 1, 
                          delta.village0 = rnorm(n.village, sd = 5),
                          delta.district0 = rnorm(n.district, sd = 5),
                          delta.cov10 = rnorm(ncol(indivcov), sd = 5),
                          lambda.cov10 = rbind(rep(0, ncol(indivcov)),
                            matrix( rnorm((K-1)*ncol(indivcov), sd = 5),
                                   nrow = K - 1, ncol = ncol(indivcov))),
                          delta.village.cov10 = rnorm(ncol(villcov), sd = 5),
                          delta.district.cov10 = rnorm(ncol(distcov), sd = 5),
                          delta.province = rnorm(n.province, sd = 5),
                          lambda.village0 = matrix(c(rep(0,n.village),
                            rep(0,n.village*(K-1))), n.village, K),
                          theta.village.cov10 = rbind(rep(0, ncol(villcov)),
                            matrix( rnorm((K-1)*ncol(villcov), sd = 5),
                                   nrow = K - 1, ncol = ncol(villcov))),
                          theta.district.cov10 = rbind(rep(0, ncol(distcov)),
                            matrix( rnorm((K-1)*ncol(distcov), sd = 5),
                                   nrow = K - 1, ncol = ncol(distcov))),
                          theta.district0 = matrix(c(rep(0,n.district),
                            rnorm(n.district*(K-1), sd = 5)), n.district, K),
                          theta.province0 = matrix(c(rep(0,n.province),
                            rnorm(n.province*(K-1), sd = 5)), n.province, K),
                          omega = c(0,runif(K-1)),
                          sigma = runif(n.district),
                          psi.district = matrix(c(rep(0,n.district),
                            runif(n.district*(K-1))), n.district, K),
                          psi.province = matrix(c(rep(0,n.province),
                            runif(n.province*(K-1))), n.province, K)
  )}
params <- c("sd.x", "beta", "alpha", "psi.district", "psi.province", "omega", "sigma.district", "sigma.province", "cor.sx", "mean.correct", "delta.village", "delta.district", "delta.province", "lambda.village", "theta.district", "theta.province", "delta.cov1", "lambda.cov1", "delta.village.cov1", "delta.district.cov1", "theta.village.cov1", "theta.district.cov1", "s") 

model <- "ModelVillageDistrictProvinceCovariates.R"

fit <- jags(data, inits, params, model.file=model, n.iter=5000, n.chains = 1)

save(fit, file = paste("../data/endorseResultsCovariates.3Aug12.",violence.exp, ".", i, ".RData", sep =""))

sink(paste("results/endorseResultsCovariates.",violence.exp, ".", i, ".ROut", sep =""))
print(fit)
sink()

check <- any(fit$BUGSoutput$summary[, "Rhat"] > 1.2) ## Set standard for convergence
counter <- 0
while(check & counter<2000){ 
	fit <- update(fit, n.iter=5000, n.thin=20)
	check <- any(fit$BUGSoutput$summary[, "Rhat"] > 1.2) ## True if any parameters have not reached convergence
	counter <- counter + 1
        
        sink(paste("results/endorseResultsCovariates.",violence.exp, ".", i, ".ROut", sep =""))
        print(fit)
        sink()

        save(fit, file = paste("results/endorseResultsCovariates.",violence.exp, ".", i, ".RData", sep =""))

}  

  
